NMDS on understory vegetation in Coastal Barrens group (MA & LI)

Libraries:

# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
library(ggordiplots)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------

Source, tranform, select, and relativize MA & LI veg data (hidden)

Species found in at least 5% (11 plots) of plots in Coastal Barrens are: VAPA, GABA, QUIL, GAPR, VAAN, PTAQ, QUCO, PIRI, KAAN, MELI, RUSP, CAPE, QUAL, SMGL, & QUPR

Source and transform envr & categorical data; merge for predictor df to compare nmds against (hidden)

Run NMDS

# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
                   distance = 'bray',
                   k = 3,
                   trymax = 20,
                   autotransform = FALSE)

## Run 0 stress 0.1272962 
## Run 1 stress 0.1275737 
## ... Procrustes: rmse 0.03524096  max resid 0.1289734 
## Run 2 stress 0.1275741 
## ... Procrustes: rmse 0.03471258  max resid 0.1272658 
## Run 3 stress 0.127296 
## ... New best solution
## ... Procrustes: rmse 0.0001272258  max resid 0.000417742 
## ... Similar to previous best
## Run 4 stress 0.1582628 
## Run 5 stress 0.1272958 
## ... New best solution
## ... Procrustes: rmse 0.000587606  max resid 0.001999951 
## ... Similar to previous best
## Run 6 stress 0.1275741 
## ... Procrustes: rmse 0.03499637  max resid 0.1277436 
## Run 7 stress 0.1272965 
## ... Procrustes: rmse 0.0008573844  max resid 0.002921016 
## ... Similar to previous best
## Run 8 stress 0.1275736 
## ... Procrustes: rmse 0.03451709  max resid 0.1262122 
## Run 9 stress 0.1275735 
## ... Procrustes: rmse 0.03479077  max resid 0.127088 
## Run 10 stress 0.1272959 
## ... Procrustes: rmse 0.0004557986  max resid 0.001556358 
## ... Similar to previous best
## Run 11 stress 0.1272965 
## ... Procrustes: rmse 0.0008645533  max resid 0.002945962 
## ... Similar to previous best
## Run 12 stress 0.1275742 
## ... Procrustes: rmse 0.03502329  max resid 0.1277945 
## Run 13 stress 0.1272961 
## ... Procrustes: rmse 0.0002420751  max resid 0.0007293543 
## ... Similar to previous best
## Run 14 stress 0.1275734 
## ... Procrustes: rmse 0.03460257  max resid 0.1264966 
## Run 15 stress 0.127574 
## ... Procrustes: rmse 0.03434109  max resid 0.1256504 
## Run 16 stress 0.1272962 
## ... Procrustes: rmse 0.0007275598  max resid 0.002468901 
## ... Similar to previous best
## Run 17 stress 0.1275736 
## ... Procrustes: rmse 0.03484616  max resid 0.1272575 
## Run 18 stress 0.1272959 
## ... Procrustes: rmse 8.512657e-05  max resid 0.000239125 
## ... Similar to previous best
## Run 19 stress 0.1272963 
## ... Procrustes: rmse 0.0007570094  max resid 0.002579365 
## ... Similar to previous best
## Run 20 stress 0.1275734 
## ... Procrustes: rmse 0.03459637  max resid 0.1264763 
## *** Best solution repeated 8 times

# converges, stress of 0.127 - 0.128

set.seed(2)
# 2 dimensions
mali_2d <- metaMDS(rel_mali5, 
                   distance = 'bray', 
                   k = 2,
                   try = 40,
                   trymax = 40,
                   autotransform = FALSE)

## Run 0 stress 0.2013098 
## Run 1 stress 0.2421891 
## Run 2 stress 0.2164813 
## Run 3 stress 0.2248403 
## Run 4 stress 0.2013098 
## ... New best solution
## ... Procrustes: rmse 7.336217e-05  max resid 0.0002404534 
## ... Similar to previous best
## Run 5 stress 0.2166909 
## Run 6 stress 0.2013098 
## ... Procrustes: rmse 2.20954e-05  max resid 7.187511e-05 
## ... Similar to previous best
## Run 7 stress 0.2299679 
## Run 8 stress 0.2334747 
## Run 9 stress 0.2067192 
## Run 10 stress 0.2063475 
## Run 11 stress 0.2248403 
## Run 12 stress 0.2052525 
## Run 13 stress 0.2013098 
## ... Procrustes: rmse 8.623799e-05  max resid 0.0002840637 
## ... Similar to previous best
## Run 14 stress 0.2248403 
## Run 15 stress 0.2013098 
## ... Procrustes: rmse 7.687422e-05  max resid 0.0002349284 
## ... Similar to previous best
## Run 16 stress 0.2385195 
## Run 17 stress 0.2378598 
## Run 18 stress 0.2013098 
## ... Procrustes: rmse 7.160904e-05  max resid 0.0002354649 
## ... Similar to previous best
## Run 19 stress 0.240332 
## Run 20 stress 0.271914 
## Run 21 stress 0.2063472 
## Run 22 stress 0.2013098 
## ... Procrustes: rmse 3.826719e-05  max resid 9.76465e-05 
## ... Similar to previous best
## Run 23 stress 0.2337596 
## Run 24 stress 0.2518641 
## Run 25 stress 0.2401184 
## Run 26 stress 0.2013098 
## ... Procrustes: rmse 6.309312e-06  max resid 1.438008e-05 
## ... Similar to previous best
## Run 27 stress 0.2431912 
## Run 28 stress 0.2013098 
## ... Procrustes: rmse 2.79146e-06  max resid 8.73469e-06 
## ... Similar to previous best
## Run 29 stress 0.2076391 
## Run 30 stress 0.2013099 
## ... Procrustes: rmse 0.0001601721  max resid 0.0005288505 
## ... Similar to previous best
## Run 31 stress 0.2013098 
## ... Procrustes: rmse 2.72951e-05  max resid 6.009303e-05 
## ... Similar to previous best
## Run 32 stress 0.2185542 
## Run 33 stress 0.2248403 
## Run 34 stress 0.2248403 
## Run 35 stress 0.239287 
## Run 36 stress 0.2248403 
## Run 37 stress 0.2076443 
## Run 38 stress 0.2394146 
## Run 39 stress 0.2013098 
## ... Procrustes: rmse 1.065455e-05  max resid 3.451836e-05 
## ... Similar to previous best
## Run 40 stress 0.2013098 
## ... Procrustes: rmse 2.107096e-05  max resid 6.965666e-05 
## ... Similar to previous best
## *** Best solution repeated 12 times

# stress still too high to represent in 2 dimensions, lowest is ~ 0.2

rm(mali_2d)

Make correlation data for bi-plots

ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)

# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))

ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)

# species data   -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)

sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)

print(sig.ma_spp.scrs)

##            NMDS1       NMDS2 Species  pval
## GAPR -0.43281215 -0.36593981    GAPR 0.028
## VAAN -0.58546249  0.23635680    VAAN 0.003
## QUCO  0.47196034  0.31546257    QUCO 0.009
## PIRI -0.09465907  0.63991646    PIRI 0.001
## KAAN -0.64342045  0.04947569    KAAN 0.001
## MELI -0.21624713 -0.54654802    MELI 0.003
## RUSP -0.62454425  0.06726456    RUSP 0.004
## SMGL  0.06314658  0.74452749    SMGL 0.001
## QUPR -0.58228890 -0.08839473    QUPR 0.012

# create species vectors for indicator species -------------------
mali.ind.spp <- ma.spp.scrs %>% 
  filter(Species %in% c('QUPR', "KAAN", "RUSP", "SMGL"))

# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)

print(ma.sig_envr.scrs)

##                 NMDS1      NMDS2 env.variables  pval
## Moss       -0.1485843  0.5469886          Moss 0.011
## BA_HA       0.5191909 -0.3902654         BA_HA 0.004
## PIRI.BA_HA  0.5939896 -0.4423759    PIRI.BA_HA 0.001

ma.sig_envr.scrs <- ma.sig_envr.scrs %>% 
  filter(env.variables != "BA_HA")

Significant species are: GAPR, VAAN, QUCO, PIRI, KAAN, MELI, RUSP, SMGL, QUPR

Significant environmental factors are: Moss, Basal area per hectare, and PIRI basal area per hectare.

I’m only going to keep one BA measure; PIRI has a lower p value, so thinking that one, but open to other ideas.

This is a 3D ordination being represented by 2D

Using the first two axes, which represent the most variation in the ordination space

library(plotly)

fort.mali3 <- fortify(mali_3d) %>% 
  filter(score == "sites")

plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, )

# Here is a 3D plot, just to show

stressplot(mali_3d)

ordicloud(mali_3d)

GG plot

# based on the idea that NMDS scores express the highest variation in first axis and on down, I will plot the 1st two and take a look -------------------
# control = yellow, fallrx = dark green, harvest = beige, mowrx = light green, spring rx = red

colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")

# create ellipse plot -------------------
mali.plot.e <- gg_ordiplot(mali_3d, groups = ma.meta.data_veg$Treat_Type)

# pull out ellipse data -------------------
ellipse.mali <- fortify(mali.plot.e$df_ellipse) %>% 
  group_by(Group)

# base plot -------------------
ma.nmds.plot <-  ggplot(ma.site.scrs, aes(x=NMDS1, y=NMDS2))+
  geom_point(aes(NMDS1, NMDS2, colour = factor(ma.site.scrs$Treat_Type), shape = factor(ma.site.scrs$Region)), size = 2)+
  scale_color_manual(values =  c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  coord_fixed()+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type", shape = "Region")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10))

# plot with ellipse -------------------
ma.nmds.plot + 
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  labs(title = "Coastal barrens ordination \n (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

# add significant species -------------------
ma.nmds.plot +
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  geom_segment(data = sig.ma_spp.scrs, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd = 0.3)+
  ggrepel::geom_text_repel(data = sig.ma_spp.scrs, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
  labs(title = "Coastal barrens ordination with \n significant species (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

# add indicators species -------------------
ind.sp.mali_color <- c("QUPR" = "#81A88D", "KAAN" = "#81A88D", "SMGL" = "#A2A475", "RUSP" = "#81A88D")

ma.nmds.plot +
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
  ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
  labs(title = "Coastal barrens ordination with \n indicator species (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

# add significant envr factors -------------------

ma.nmds.plot +
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  geom_segment(data = ma.sig_envr.scrs, aes(x = 0, xend = NMDS1, y = 0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")))+
  ggrepel::geom_text_repel(data = ma.sig_envr.scrs, aes(x=NMDS1, y=NMDS2, label = env.variables), cex = 4, direction = "both", segment.size = 0.25) +
  labs(title = "Coastal barrens ordination with \n environmental vectors (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

Plot with significant species & envr factors

ma.nmds.plot +
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  geom_segment(data = sig.ma_spp.scrs, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd = 0.3)+
  ggrepel::geom_text_repel(data = sig.ma_spp.scrs, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
  geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
  ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
  labs(title = "Coastal barrens ordination with significant \n species & environmental factors (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

Plot with indcator species & envr factors

ma.nmds.plot +
  geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group),  linejoin = "round", show.legend = FALSE)+
  geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
  ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
  geom_segment(data = ma.sig_envr.scrs, aes(x = 0, xend = NMDS1, y = 0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")))+
  ggrepel::geom_text_repel(data = ma.sig_envr.scrs, aes(x=NMDS1, y=NMDS2, label = env.variables), cex = 4, direction = "both", segment.size = 0.25) +
  labs(title = "Coastal barrens ordination with indicator \n species & environmental factors (Axes 1 & 2)")+
  theme(plot.title=element_text(hjust=0.5, size=20))

Run PerMANOVA to look at the impacts of region and treatment on community composition

summary(as.factor(ma.meta.data_veg$Treat_Type))

##  Control   FallRx  Harvest    MowRx SpringRx 
##        5        3        6        3        6

ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)

print(ma.tt)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
##            Df SumOfSqs      R2      F Pr(>F)    
## Treat_Type  4   1.9569 0.31903 2.1082  0.001 ***
## Residual   18   4.1769 0.68097                  
## Total      22   6.1339 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise follow up

library(pairwiseAdonis)

pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)

## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
## 
## $FallRx_vs_SpringRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.21549 0.10956 0.8613  0.614
## Residual    7  1.75138 0.89044              
## Total       8  1.96687 1.00000              
## 
## $FallRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.48859 0.46035 3.4122    0.1
## Residual    4  0.57276 0.53965              
## Total       5  1.06134 1.00000              
## 
## $FallRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.40379 0.20243 1.7767  0.109
## Residual    7  1.59094 0.79757              
## Total       8  1.99473 1.00000              
## 
## $FallRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.39033 0.23712 1.8649   0.12
## Residual    6  1.25581 0.76288              
## Total       7  1.64614 1.00000              
## 
## $SpringRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.65697 0.28336 2.7678  0.015 *
## Residual    7  1.66151 0.71664                
## Total       8  2.31848 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1   0.4952 0.15598 1.8481  0.061 .
## Residual   10   2.6797 0.84402                
## Total      11   3.1749 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Control
##            Df SumOfSqs     R2      F Pr(>F)
## Treat_Type  1  0.16729 0.0666 0.6422  0.778
## Residual    9  2.34457 0.9334              
## Total      10  2.51186 1.0000              
## 
## $MowRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.72529 0.32578 3.3823  0.014 *
## Residual    7  1.50107 0.67422                
## Total       8  2.22636 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $MowRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.84147 0.41918 4.3303  0.012 *
## Residual    6  1.16594 0.58082                
## Total       7  2.00741 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Harvest_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.58418 0.21102 2.4072  0.023 *
## Residual    9  2.18412 0.78898                
## Total      10  2.76830 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

#significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control

Indicator species analysis

library(labdsv)

ma.meta.data_veg$Treat_Group <- NA

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)

mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)

summary(mali_isa)

##      cluster indicator_value probability
## QUPR       2          0.8549       0.001
## KAAN       2          0.7534       0.012
## RUSP       2          0.7097       0.006
## SMGL       4          0.6896       0.009
## 
## Sum of probabilities                 =  4.689 
## 
## Sum of Indicator Values              =  6.68 
## 
## Sum of Significant Indicator Values  =  3.01 
## 
## Number of Significant Indicators     =  4 
## 
## Significant Indicator Distribution
## 
## 2 4 
## 3 1

gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]

print(fridg)

##      group    indval pvalue freq
## QUPR     2 0.8549332  0.001    6
## KAAN     2 0.7534400  0.012   10
## RUSP     2 0.7096551  0.006    7
## SMGL     4 0.6896051  0.009    6

# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is Control